In [6]:
aux_plt = az.plot_khat(res, show_bins=True, figsize=(7,5))
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mplt
import seaborn as sns
import arviz as az
import datetime
import fit_arima
import open_data
df = open_data.open()
ritorno = fit_arima.compute(2, 1, catene=4, samples_per_chain=2500, burnin=1000)
21:10:52 - cmdstanpy - INFO - CmdStan start processing
chain 1 | | 00:00 Status
chain 2 | | 00:00 Status
chain 3 | | 00:00 Status
chain 4 | | 00:00 Status
21:38:05 - cmdstanpy - INFO - CmdStan done processing. 21:38:05 - cmdstanpy - WARNING - Non-fatal error during sampling: Exception: normal_lpdf: Scale parameter is -0.111703, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61) Exception: code_model_namespace::log_prob: phi[2][11] is 1, but must be less than or equal to 1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 83, column 2 to column 47) Exception: code_model_namespace::log_prob: phi[1][1] is -nan, but must be greater than or equal to -1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 83, column 2 to column 47) Exception: code_model_namespace::log_prob: phi[1][33] is 1, but must be less than or equal to 1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 83, column 2 to column 47) Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 157, column 4 to column 34) Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 157, column 4 to column 34) Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 157, column 4 to column 34) Exception: code_model_namespace::log_prob: phi[2][8] is -nan, but must be greater than or equal to -1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 83, column 2 to column 47) Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 157, column 4 to column 34) Exception: normal_lpdf: Scale parameter is -0.00468339, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61) Exception: normal_lpdf: Scale parameter is -0.00781159, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61) Exception: normal_lpdf: Scale parameter is -19.3346, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61) Exception: normal_lpdf: Scale parameter is -20.4458, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61) Exception: normal_lpdf: Scale parameter is -1.84055, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61) Exception: code_model_namespace::log_prob: phi[2][4] is -nan, but must be greater than or equal to -1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 83, column 2 to column 47) Exception: normal_lpdf: Scale parameter is -9.85741, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61) Exception: normal_lpdf: Scale parameter is -0.00793045, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61) Exception: normal_lpdf: Scale parameter is -1.48033, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61) Exception: normal_lpdf: Scale parameter is -0.370182, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61) Exception: normal_lpdf: Scale parameter is -3.13815, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61) Exception: normal_lpdf: Scale parameter is -2.7382, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61) Exception: normal_lpdf: Scale parameter is -0.274406, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61) Exception: code_model_namespace::log_prob: phi[1][10] is -nan, but must be greater than or equal to -1.000000 (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 83, column 2 to column 47) Exception: normal_lpdf: Scale parameter is -3.44707, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61) Exception: normal_lpdf: Scale parameter is -0.00909615, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61) Exception: normal_lpdf: Scale parameter is -0.0018971, but must be positive! (in '/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/code.stan', line 143, column 4 to column 61) Consider re-running with show_console=True if the above output is unclear!
/home/br1/PythonProjects/PythonStats/ARIMA_last_computation/NoCovariates/fit_arima.py:37: FutureWarning: Indexing with a float is deprecated, and will raise an IndexError in pandas 2.0. You can manually convert to an integer key instead.
y_missing_list.append({'Stazione':df.columns[v[i]],'Data':df.index[u[i]],'Samples':inference_data.posterior.w.values[:,:,i].reshape(catene*samples_per_chain)})
Tempo di computazione: 27:12
res = az.loo(ritorno['inference_data'], var_name="log_lik", pointwise=True, scale='deviance')
res
/home/br1/PythonProjects/PythonStats/.venv/lib/python3.10/site-packages/arviz/stats/stats.py:803: UserWarning: Estimated shape parameter of Pareto distribution is greater than 0.7 for one or more samples. You should consider using a more robust model, this is because importance sampling is less likely to work well if the marginal posterior and LOO posterior are very different. This is more likely to happen with a non-robust model and highly influential observations. warnings.warn(
Computed from 10000 posterior samples and 17191 observations log-likelihood matrix.
Estimate SE
deviance_loo -14506.23 309.87
p_loo 296.12 -
There has been a warning during the calculation. Please check the results.
------
Pareto k diagnostic values:
Count Pct.
(-Inf, 0.5] (good) 17112 99.5%
(0.5, 0.7] (ok) 68 0.4%
(0.7, 1] (bad) 11 0.1%
(1, Inf) (very bad) 0 0.0%
aux_plt = az.plot_khat(res, show_bins=True, figsize=(7,5))
res = az.waic(ritorno['inference_data'], var_name="log_lik", scale='deviance')
res
/home/br1/PythonProjects/PythonStats/.venv/lib/python3.10/site-packages/arviz/stats/stats.py:1645: UserWarning: For one or more samples the posterior variance of the log predictive densities exceeds 0.4. This could be indication of WAIC starting to fail. See http://arxiv.org/abs/1507.04544 for details warnings.warn(
Computed from 10000 posterior samples and 17191 observations log-likelihood matrix.
Estimate SE
deviance_waic -14528.97 309.61
p_waic 284.76 -
There has been a warning during the calculation. Please check the results.
plt_aux = az.plot_trace(ritorno['inference_data'], var_names=['phi','theta','hyper_gamma_phi_m','hyper_gamma_phi_s'], divergences=True, figsize=(15,15))
plt_aux = az.plot_trace(ritorno['inference_data'], var_names=['mh_sigma','sh_sigma','sigma'], divergences=True, figsize=(15,10))
plt_aux = az.plot_trace(ritorno['inference_data'], var_names=['hyper_y_start_m','hyper_y_start_s','c'], divergences=True, figsize=(15,10))
az.summary(ritorno['inference_data'],var_names=['y_start'])
| mean | sd | hdi_3% | hdi_97% | mcse_mean | mcse_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|
| y_start[0, 0] | 1.422 | 0.089 | 1.251 | 1.587 | 0.002 | 0.001 | 2935.0 | 3563.0 | 1.0 |
| y_start[0, 1] | 1.419 | 0.091 | 1.248 | 1.591 | 0.002 | 0.001 | 3125.0 | 3114.0 | 1.0 |
| y_start[0, 2] | 1.422 | 0.090 | 1.252 | 1.593 | 0.002 | 0.001 | 3105.0 | 4272.0 | 1.0 |
| y_start[0, 3] | 1.424 | 0.090 | 1.253 | 1.594 | 0.002 | 0.001 | 3232.0 | 3953.0 | 1.0 |
| y_start[0, 4] | 1.427 | 0.092 | 1.245 | 1.595 | 0.002 | 0.001 | 2847.0 | 3181.0 | 1.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| y_start[2, 44] | 1.445 | 0.054 | 1.345 | 1.548 | 0.001 | 0.001 | 1875.0 | 2615.0 | 1.0 |
| y_start[2, 45] | 1.442 | 0.055 | 1.340 | 1.549 | 0.001 | 0.001 | 1967.0 | 2244.0 | 1.0 |
| y_start[2, 46] | 1.438 | 0.054 | 1.334 | 1.535 | 0.001 | 0.001 | 2333.0 | 3204.0 | 1.0 |
| y_start[2, 47] | 1.452 | 0.055 | 1.348 | 1.556 | 0.001 | 0.001 | 1877.0 | 2748.0 | 1.0 |
| y_start[2, 48] | 1.459 | 0.057 | 1.360 | 1.572 | 0.001 | 0.001 | 1908.0 | 2766.0 | 1.0 |
147 rows × 9 columns
az.summary(ritorno['inference_data'], var_names=['hyper_y_start_m', 'hyper_y_start_s'])
| mean | sd | hdi_3% | hdi_97% | mcse_mean | mcse_sd | ess_bulk | ess_tail | r_hat | |
|---|---|---|---|---|---|---|---|---|---|
| hyper_y_start_m | 1.432 | 0.038 | 1.364 | 1.506 | 0.002 | 0.001 | 578.0 | 852.0 | 1.00 |
| hyper_y_start_s | 0.082 | 0.012 | 0.061 | 0.106 | 0.000 | 0.000 | 688.0 | 1692.0 | 1.01 |
aux_shape = ritorno['inference_data'].posterior.gamma_th.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.gamma_th.values.reshape((aux_shape[0]*aux_shape[1])), legend=False)
<AxesSubplot: ylabel='Density'>
aux_shape = ritorno['inference_data'].posterior.phi.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.phi.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2],aux_shape[3]))[:,0,:], legend=False)
<AxesSubplot: ylabel='Density'>
aux_shape = ritorno['inference_data'].posterior.phi.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.phi.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2],aux_shape[3]))[:,1,:], legend=False)
<AxesSubplot: ylabel='Density'>
aux_shape = ritorno['inference_data'].posterior.c.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.c.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)
<AxesSubplot: ylabel='Density'>
aux_shape = ritorno['inference_data'].posterior.annual_mean.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.annual_mean.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)
<AxesSubplot: ylabel='Density'>
aux_shape = ritorno['inference_data'].posterior.annual_median.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.annual_median.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)
<AxesSubplot: ylabel='Density'>
aux_shape = ritorno['inference_data'].posterior.annual_max.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.annual_max.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)
<AxesSubplot: ylabel='Density'>
valori_31_dic = []
for dizionario in ritorno['missing']:
if dizionario['Data'] == '2018-12-31':
valori_31_dic.append(dizionario['Samples'])
res_plot = sns.kdeplot(valori_31_dic, legend=False)
aux_shape = ritorno['inference_data'].posterior.annual_days_over_threshold.values.shape
sns.kdeplot(ritorno['inference_data'].posterior.annual_days_over_threshold.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), legend=False)
<AxesSubplot: ylabel='Density'>
pd.set_option('display.max_columns', 500)
aux_results = pd.DataFrame(ritorno['inference_data'].posterior.annual_days_over_threshold.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), columns=df.columns).describe()
aux_results.sort_values('mean', axis=1, ascending=True)
| CASTELLUCCIO | CORTE BRUGNATELLA | SAVIGNANO DI RIGO | FEBBIO | SAN LEO | VERUCCHIO | BADIA | GIARDINI MARGHERITA | MARECCHIA | SAN PIETRO CAPOFIUME | PARCO BERTOZZI | DELTA CERVIA | LUGAGNANO | DE AMICIS | SAN LAZZARO | PARCO RESISTENZA | FRANCHINI-ANGELONI | GHERARDI | VIA CHIARINI | GAVELLO | SAVIGNANO | BESENZONE | ROMA | VILLA FULVIA | PORTA SAN FELICE | CAORLE | CASTELLARANO | CENTO | ZALAMELLA | PARCO MONTECUCCO | PARCO EDILCARANI | BOGOLESE | SARAGAT | PARADIGNA | MALCANTONE | S. LAZZARO | REMESINA | S. ROCCO | CENO | ISONZO | PARCO FERRARI | GIORDANI-FARNESE | CITTADELLA | SAN FRANCESCO | FLAMINIA | MONTEBELLO | GIARDINI | GERBIDO | TIMAVO | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| count | 10000.00000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.00000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.00000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.00000 | 10000.00000 | 10000.00000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.00000 | 10000.000000 | 10000.000000 | 10000.000000 |
| mean | 0.10200 | 0.115000 | 0.118500 | 0.302400 | 1.567600 | 4.877900 | 7.790400 | 13.415600 | 15.608700 | 15.637100 | 17.239500 | 17.616200 | 17.707000 | 17.769500 | 18.701500 | 19.475400 | 21.104400 | 21.166500 | 21.361500 | 23.260100 | 23.932500 | 23.949300 | 28.18090 | 28.534400 | 28.738900 | 29.321800 | 29.62880 | 30.606900 | 30.616000 | 30.740400 | 31.929900 | 32.091500 | 34.509200 | 35.694200 | 35.87640 | 36.64660 | 36.97600 | 38.350700 | 41.990600 | 42.731600 | 42.800400 | 46.467400 | 46.507200 | 50.969500 | 51.176500 | 54.85700 | 60.543300 | 63.666400 | 70.191400 |
| std | 0.36443 | 0.394703 | 0.398087 | 0.680147 | 1.962961 | 4.054311 | 5.462239 | 8.325383 | 9.067712 | 9.744521 | 10.226974 | 9.588717 | 9.695842 | 9.423346 | 10.193998 | 10.424829 | 11.205492 | 11.605751 | 11.138134 | 11.271324 | 11.896567 | 12.455571 | 13.45773 | 13.180083 | 13.040335 | 14.381499 | 12.01704 | 14.105356 | 14.220372 | 14.289571 | 12.852171 | 13.487006 | 14.189815 | 15.221756 | 15.03205 | 15.35704 | 14.94129 | 15.436954 | 16.984556 | 16.497745 | 16.676348 | 17.410958 | 17.140721 | 17.009283 | 18.785703 | 18.80313 | 19.254184 | 22.584466 | 21.885408 |
| min | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 1.00000 | 0.000000 | 0.000000 | 0.000000 | 2.00000 | 1.000000 | 1.000000 | 0.000000 | 2.000000 | 1.000000 | 1.000000 | 2.000000 | 0.00000 | 2.00000 | 2.00000 | 3.000000 | 4.000000 | 1.000000 | 1.000000 | 4.000000 | 7.000000 | 9.000000 | 8.000000 | 8.00000 | 10.000000 | 5.000000 | 9.000000 |
| 25% | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 2.000000 | 4.000000 | 7.000000 | 9.000000 | 9.000000 | 10.000000 | 11.000000 | 11.000000 | 11.000000 | 11.000000 | 12.000000 | 13.000000 | 13.000000 | 13.000000 | 15.000000 | 15.000000 | 15.000000 | 19.00000 | 19.000000 | 19.000000 | 19.000000 | 21.00000 | 20.000000 | 20.000000 | 20.000000 | 23.000000 | 22.000000 | 24.000000 | 25.000000 | 25.00000 | 26.00000 | 26.00000 | 27.000000 | 30.000000 | 31.000000 | 31.000000 | 34.000000 | 34.000000 | 39.000000 | 38.000000 | 41.00000 | 47.000000 | 48.000000 | 55.000000 |
| 50% | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 1.000000 | 4.000000 | 7.000000 | 12.000000 | 14.000000 | 14.000000 | 15.000000 | 16.000000 | 16.000000 | 16.000000 | 17.000000 | 18.000000 | 19.000000 | 19.000000 | 20.000000 | 21.000000 | 22.000000 | 22.000000 | 26.00000 | 27.000000 | 27.000000 | 27.000000 | 28.00000 | 29.000000 | 28.000000 | 29.000000 | 30.000000 | 31.000000 | 33.000000 | 34.000000 | 34.00000 | 35.00000 | 35.00000 | 36.000000 | 40.000000 | 41.000000 | 41.000000 | 45.000000 | 45.000000 | 50.000000 | 49.000000 | 53.00000 | 59.000000 | 61.000000 | 69.000000 |
| 75% | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 2.000000 | 7.000000 | 10.000000 | 18.000000 | 20.250000 | 21.000000 | 22.000000 | 23.000000 | 23.000000 | 23.000000 | 24.000000 | 25.000000 | 27.000000 | 27.000000 | 27.000000 | 30.000000 | 31.000000 | 31.000000 | 35.00000 | 36.000000 | 36.000000 | 37.000000 | 37.00000 | 38.250000 | 39.000000 | 39.000000 | 39.000000 | 40.000000 | 43.000000 | 44.000000 | 44.25000 | 45.00000 | 46.00000 | 47.000000 | 52.000000 | 52.000000 | 53.000000 | 57.000000 | 57.000000 | 61.000000 | 62.000000 | 66.00000 | 72.000000 | 77.000000 | 83.000000 |
| max | 4.00000 | 7.000000 | 6.000000 | 8.000000 | 21.000000 | 38.000000 | 41.000000 | 80.000000 | 84.000000 | 92.000000 | 91.000000 | 83.000000 | 86.000000 | 85.000000 | 83.000000 | 104.000000 | 102.000000 | 112.000000 | 85.000000 | 101.000000 | 137.000000 | 105.000000 | 130.00000 | 166.000000 | 101.000000 | 111.000000 | 98.00000 | 121.000000 | 119.000000 | 124.000000 | 100.000000 | 119.000000 | 120.000000 | 130.000000 | 130.00000 | 144.00000 | 126.00000 | 120.000000 | 153.000000 | 156.000000 | 141.000000 | 155.000000 | 157.000000 | 139.000000 | 159.000000 | 156.00000 | 163.000000 | 199.000000 | 218.000000 |
pd.set_option('display.max_columns', 500)
aux_results = pd.DataFrame(ritorno['inference_data'].posterior.is_over_daily_limit.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), columns=df.columns).describe()
aux_results.sort_values('mean', axis=1, ascending=True)
| CASTELLUCCIO | FEBBIO | CORTE BRUGNATELLA | SAVIGNANO DI RIGO | SAN LEO | VERUCCHIO | BADIA | GIARDINI MARGHERITA | MARECCHIA | SAN PIETRO CAPOFIUME | DE AMICIS | LUGAGNANO | DELTA CERVIA | PARCO BERTOZZI | SAN LAZZARO | PARCO RESISTENZA | FRANCHINI-ANGELONI | GHERARDI | VIA CHIARINI | GAVELLO | SAVIGNANO | BESENZONE | ROMA | VILLA FULVIA | PORTA SAN FELICE | CASTELLARANO | CAORLE | CENTO | PARCO MONTECUCCO | ZALAMELLA | PARCO EDILCARANI | BOGOLESE | SARAGAT | PARADIGNA | MALCANTONE | S. LAZZARO | REMESINA | S. ROCCO | CENO | ISONZO | PARCO FERRARI | GIORDANI-FARNESE | CITTADELLA | FLAMINIA | SAN FRANCESCO | MONTEBELLO | GERBIDO | GIARDINI | TIMAVO | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| count | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.00000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.00000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.00000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.00000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 |
| mean | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.000200 | 0.001300 | 0.019700 | 0.034300 | 0.041900 | 0.045500 | 0.051700 | 0.05210 | 0.057100 | 0.066200 | 0.078100 | 0.105000 | 0.108900 | 0.111100 | 0.135500 | 0.157300 | 0.157400 | 0.249100 | 0.258900 | 0.26580 | 0.275400 | 0.287700 | 0.313500 | 0.316800 | 0.319500 | 0.34190 | 0.356600 | 0.426600 | 0.452100 | 0.453500 | 0.47890 | 0.485400 | 0.524300 | 0.610400 | 0.635600 | 0.638300 | 0.717900 | 0.721100 | 0.794600 | 0.816300 | 0.858800 | 0.912300 | 0.916200 | 0.961700 |
| std | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.014141 | 0.036034 | 0.138974 | 0.182008 | 0.200371 | 0.208408 | 0.221432 | 0.22224 | 0.232045 | 0.248644 | 0.268342 | 0.306569 | 0.311529 | 0.314272 | 0.342274 | 0.364102 | 0.364196 | 0.432513 | 0.438052 | 0.44178 | 0.446738 | 0.452713 | 0.463939 | 0.465252 | 0.466306 | 0.47437 | 0.479019 | 0.494608 | 0.497725 | 0.497858 | 0.49958 | 0.499812 | 0.499434 | 0.487684 | 0.481286 | 0.480517 | 0.450044 | 0.448481 | 0.404014 | 0.387259 | 0.348245 | 0.282872 | 0.277101 | 0.191929 |
| min | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| 25% | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 |
| 50% | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.00000 | 0.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 |
| 75% | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.00000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 1.000000 | 1.00000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.00000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.00000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 |
| max | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.00000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.00000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.00000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.00000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 |
aux_results.loc['mean',:].plot(kind='barh', figsize=(15,15), xticks=np.linspace(0,1,11))
#res = plt.xticks(rotation = 30)
<AxesSubplot: >
pd.set_option('display.max_columns', 500)
aux_results = pd.DataFrame(ritorno['inference_data'].posterior.is_over_annual_limit.values.reshape((aux_shape[0]*aux_shape[1],aux_shape[2])), columns=df.columns).describe()
aux_results.sort_values('mean', axis=1, ascending=True)
| CASTELLUCCIO | ROMA | CASTELLARANO | PARCO EDILCARANI | DELTA CERVIA | BESENZONE | PORTA SAN FELICE | GAVELLO | ZALAMELLA | PARCO MONTECUCCO | CENTO | BOGOLESE | SARAGAT | REMESINA | S. ROCCO | GHERARDI | LUGAGNANO | CAORLE | SAN PIETRO CAPOFIUME | FEBBIO | CORTE BRUGNATELLA | SAN LAZZARO | SAN LEO | VERUCCHIO | BADIA | GIARDINI MARGHERITA | SAVIGNANO DI RIGO | PARCO BERTOZZI | MARECCHIA | DE AMICIS | VIA CHIARINI | FRANCHINI-ANGELONI | PARCO RESISTENZA | VILLA FULVIA | MALCANTONE | SAVIGNANO | S. LAZZARO | CENO | ISONZO | PARADIGNA | SAN FRANCESCO | PARCO FERRARI | GIORDANI-FARNESE | CITTADELLA | FLAMINIA | MONTEBELLO | GIARDINI | GERBIDO | TIMAVO | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| count | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0 | 10000.0000 | 10000.0000 | 10000.0000 | 10000.0000 | 10000.0000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 | 10000.000000 |
| mean | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0001 | 0.0001 | 0.0001 | 0.0001 | 0.0001 | 0.000200 | 0.000200 | 0.000300 | 0.000400 | 0.000400 | 0.000600 | 0.001400 | 0.001700 | 0.002400 | 0.014100 | 0.014600 |
| std | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0100 | 0.0100 | 0.0100 | 0.0100 | 0.0100 | 0.014141 | 0.014141 | 0.017319 | 0.019997 | 0.019997 | 0.024489 | 0.037392 | 0.041198 | 0.048933 | 0.117909 | 0.119951 |
| min | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| 25% | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| 50% | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| 75% | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
| max | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 | 1.000000 |
aux_results.loc['mean',:].plot(kind='barh', figsize=(15,15), xticks=np.linspace(0,1,11))
#res = plt.xticks(rotation = 30)
<AxesSubplot: >
y_post_pred = ritorno['inference_data'].posterior.y_post_pred.to_numpy()
y_post_pred = y_post_pred.reshape(y_post_pred.shape[0]*y_post_pred.shape[1], y_post_pred.shape[2], y_post_pred.shape[3])
y_post_pred.shape
(10000, 365, 49)
import ipywidgets as widgets
from ipywidgets import HBox, VBox
from IPython.display import display
@widgets.interact(stazione = df.columns)
def f(stazione):
col_map = sns.light_palette((20, 75, 70), input='husl', as_cmap=True)
plt.figure(figsize=(14,8))
ax = plt.subplot(1,1,1)
idx_stazione = df.columns.to_list().index(stazione)
"""
sns.lineplot(np.transpose(y_post_pred[0:50,:,idx_stazione]))
"""
lower_lim = np.zeros(len(df.index))
upper_lim = np.zeros(len(df.index))
for i in range(len(df.index)):
lower_lim[i] = np.percentile(y_post_pred[:,i,idx_stazione],2.5)
upper_lim[i] = np.percentile(y_post_pred[:,i,idx_stazione],97.5)
# sns.lineplot(lower_lim)
# sns.lineplot(upper_lim)
for i in range(len(df.index)):
ax.add_patch(mplt.patches.Rectangle((i,lower_lim[i]),1,upper_lim[i]-lower_lim[i], fill=True, color=col_map(180)))
sns.lineplot(df[stazione])
sns.lineplot(np.mean(y_post_pred[:,:,idx_stazione], axis=0))
index = 0
for line in ax.get_lines():
if index == 0:
col_map = sns.dark_palette((230,90,65), input='husl', as_cmap=True)
line.set_c(col_map(155))
else:
col_map = sns.dark_palette((10,90,65), input='husl', as_cmap=True)
line.set_c(col_map(175))
index += 1
plt.ylim(bottom=0,top=2.5)
primo_giorno = datetime.date(2018,1,1)
date_da_segnare = []
date_da_segnare_posizioni = []
for i in range(12):
date_da_segnare.append(datetime.date(2018,i+1,1))
date_da_segnare_posizioni.append((date_da_segnare[2*i] - primo_giorno).days)
date_da_segnare.append(datetime.date(2018,i+1,15))
date_da_segnare_posizioni.append((date_da_segnare[2*i+1] - primo_giorno).days)
date_da_segnare[2*i] = date_da_segnare[2*i].isoformat()
date_da_segnare[2*i+1] = date_da_segnare[2*i+1].isoformat()
plt.xticks(date_da_segnare_posizioni,date_da_segnare,rotation=65)
plt.tick_params(
axis='x', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
bottom=True, # ticks along the bottom edge are off
top=False, # ticks along the top edge are off
labelbottom=True)
plt.grid()
"""
ax.get_legend().remove()
col_vals = np.linspace(1,255,num=len(df.columns))
index = 0
for line in ax.get_lines():
if(index == len(ax.get_lines()) - 1):
col_map = sns.dark_palette((230,90,65), input='husl', as_cmap=True)
line.set_c(col_map(175))
continue
if(index == len(ax.get_lines()) - 2):
col_map = sns.dark_palette((120,90,65), input='husl', as_cmap=True)
line.set_c(col_map(175))
index += 1
continue
#line.set_c(col_map(int(np.round(col_vals[index]))))
line.set_c(col_map(220))
line.set_alpha(0.3)
index += 1
"""
plt.show()
interactive(children=(Dropdown(description='stazione', options=('CASTELLUCCIO', 'FEBBIO', 'CORTE BRUGNATELLA',…
def f(stazione):
col_map = sns.light_palette((20, 75, 70), input='husl', as_cmap=True)
plt.figure(figsize=(14,8))
ax = plt.subplot(1,1,1)
idx_stazione = df.columns.to_list().index(stazione)
"""
sns.lineplot(np.transpose(y_post_pred[0:50,:,idx_stazione]))
"""
lower_lim = np.zeros(len(df.index))
upper_lim = np.zeros(len(df.index))
for i in range(len(df.index)):
lower_lim[i] = np.percentile(y_post_pred[:,i,idx_stazione],2.5)
upper_lim[i] = np.percentile(y_post_pred[:,i,idx_stazione],97.5)
# sns.lineplot(lower_lim)
# sns.lineplot(upper_lim)
for i in range(len(df.index)):
ax.add_patch(mplt.patches.Rectangle((i,lower_lim[i]),1,upper_lim[i]-lower_lim[i], fill=True, color=col_map(180)))
sns.lineplot(df[stazione])
sns.lineplot(np.mean(y_post_pred[:,:,idx_stazione], axis=0))
index = 0
for line in ax.get_lines():
if index == 0:
col_map = sns.dark_palette((230,90,65), input='husl', as_cmap=True)
line.set_c(col_map(155))
else:
col_map = sns.dark_palette((10,90,65), input='husl', as_cmap=True)
line.set_c(col_map(175))
index += 1
plt.ylim(bottom=0,top=2.5)
primo_giorno = datetime.date(2018,1,1)
date_da_segnare = []
date_da_segnare_posizioni = []
for i in range(12):
date_da_segnare.append(datetime.date(2018,i+1,1))
date_da_segnare_posizioni.append((date_da_segnare[2*i] - primo_giorno).days)
date_da_segnare.append(datetime.date(2018,i+1,15))
date_da_segnare_posizioni.append((date_da_segnare[2*i+1] - primo_giorno).days)
date_da_segnare[2*i] = date_da_segnare[2*i].isoformat()
date_da_segnare[2*i+1] = date_da_segnare[2*i+1].isoformat()
plt.xticks(date_da_segnare_posizioni,date_da_segnare,rotation=65)
plt.tick_params(
axis='x', # changes apply to the x-axis
which='both', # both major and minor ticks are affected
bottom=True, # ticks along the bottom edge are off
top=False, # ticks along the top edge are off
labelbottom=True)
plt.grid()
"""
ax.get_legend().remove()
col_vals = np.linspace(1,255,num=len(df.columns))
index = 0
for line in ax.get_lines():
if(index == len(ax.get_lines()) - 1):
col_map = sns.dark_palette((230,90,65), input='husl', as_cmap=True)
line.set_c(col_map(175))
continue
if(index == len(ax.get_lines()) - 2):
col_map = sns.dark_palette((120,90,65), input='husl', as_cmap=True)
line.set_c(col_map(175))
index += 1
continue
#line.set_c(col_map(int(np.round(col_vals[index]))))
line.set_c(col_map(220))
line.set_alpha(0.3)
index += 1
"""
plt.show()
for stazione in df.columns:
f(stazione)
